#ifndef _EQUATION_SOLVER_H_
#define _EQUATION_SOLVER_H_
#include <cmath>
#include <iostream>
#include <limits>
#include <algorithm>

const double eps=1e-8;


class Function{
    public:
	    virtual double operator()(const double &x)  = 0;
        virtual double diff ( const double &x) {
            return ((*this)(x+eps) - (*this)(x-eps)) / (2*eps);
    }
        virtual double diff2 ( const double &x) {
            return ((*this)(x+eps) -(*this)(x)+ (*this)(x-eps)) / (eps*eps);
    }
};


class EquationSolver{
    protected:
        virtual double solve()=0;
};

class BisectionSolver :public EquationSolver {
    private:
        Function &f;
        double a,b,delta;
        int M;
    public:
        BisectionSolver(Function &f,double a,double b, double delta, int M ): f(f), a(a), b(b), delta(delta), M(M) {};
        virtual double solve(){
            int k;
            double h = b - a;
            double u = f(a);
            double w,c;
            for(k = 0; k <= M; k++){
                h=h/2;
                c = a + h ;
                w = f(c);
                if( fabs(h) < delta || fabs(w) < eps )
                    break;
                else if( w > 0 && u > 0 || w < 0 && u < 0 )
                    a=c;
            }
            std::cerr << "Bisection" << std::endl;
            std::cerr << "Total Step : "<< k << std::endl;
            std::cerr << "Approxiate Root: " << c << std::endl;
            std::cerr << "Function Value: " << f(c) << std::endl;
            return c;
        }
};


class NewtonSolver : public EquationSolver{
    private:
        Function &f;
        double x0;
        int M;
    public:
        NewtonSolver(Function &f,double x0,int M):f(f), x0(x0),  M(M) {}
        virtual double solve(){
            int k;
            double x=x0;
            double u;
            for(k = 0; k <= M; k++){
                u=f(x);
                if( fabs(u) < eps )
                    break;
                else 
                    x=x-u/f.diff(x);
            }
            std::cerr <<"Newton" << std::endl;
            std::cerr <<"Total Step : "<< k << std::endl;
            std::cerr << "Initial Point: " << x0 << std::endl;
            std::cerr <<"Approxiate Root: " << x << std::endl;
            std::cerr << "Function Value: " << f(x) << std::endl;
            return x;
        }
};

class SecantSolver:public EquationSolver{
    private:
        Function &f;
        double a,b,delta;
        int M;
    public:
        SecantSolver(Function &f,double a,double b,double delta,int M):f(f), a(a), b(b), delta(delta), M(M) {};
        virtual double solve(){
            double x0=a,x1=b;
            double u=f(x1);
            double v=f(x0);
            double s;
            int k=2;
            for (k = 2; k <= M; k++){
                if (fabs(u) > fabs(v)) {
                    std::swap(x0, x1);
                    std::swap(u, v);
                }
                s = (x1 - x0) / (u - v);
                x0 = x1;
                v = u;
                x1 = x1 - u * s;
                u = f(x1);
                u = f(x1);
                if( fabs(x1 - x0) < delta || fabs(u) < eps )
                    break;
            }
            std::cerr <<"Secant" << std::endl;
            std::cerr <<"Total Step : "<< k << std::endl;
            std::cerr <<"Approxiate Root: " << x1 << std::endl;
            std::cerr << "Function Value: " << f(x1) << std::endl;
            return x1;
        }
};
#endif